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Abstract 



Models of inflationary cosmology can lead to variation of observable pa- 
rameters ("constants of Nature") on extremely large scales. The question of 
making probabilistic predictions for today's observables in such models has 
been investigated in the literature. Because of the infinite thermalized volume 
resulting from eternal inflation, it has proven difficult to obtain a meaningful 
and unambiguous probability distribution for observables, in particular due 
to the gauge dependence. In the present paper, we further develop the gauge- 
invariant procedure proposed in a previous work for models with a continuous 
variation of "constants" . The recipe uses an unbiased selection of a connected 
piece of the thermalized volume as sample for the probability distribution. To 
implement the procedure numerically, we develop two methods applicable to 
a reasonably wide class of models: one based on the Fokker-Planck equation 
of stochastic inflation, and the other based on direct simulation of inflationary 
spacetime. We present and compare results obtained using these methods. 



I. INTRODUCTION 

The parameters we call "constants of Nature" may in fact be variables related to some 
slowly varying fields. For example, what we perceive as a cosmological constant could be a 
potential U(x) of some field x( x )- If this potential is very flat, so that the evolution of x 
is much slower than the Hubble expansion, then observations will not distinguish between 
U(x) an d a true cosmological constant. Observers in different parts of the universe could 
then measure different values of U(x)- 

Spatial variation of the fields Xa associated with the "constants" can naturally arise 
in the framework of inflationary cosmology [|TJ. The dynamics of light scalar fields during 
inflation are strongly influenced by quantum fluctuations, so different regions of the universe 
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thermalize with different values of Xa- An intriguing question is whether or not we can predict 
the values of the "constants" we are most likely to observe. 

It is important to realize that this question arises not only in some exotic models es- 
pecially designed to have variable constants. On the contrary, we have to face it in all 
inflationary models attempting to predict the spectrum of cosmological density fluctuations. 
The density fluctuation 5p/p(l) is determined by the quantum fluctuation 5(f)(1) of the infla- 
ton field <p at the time when the corresponding length scale I crossed the horizon. Different 
realizations of random fluctuations 6(f)(1) result in different density fluctuation spectra in 
widely separated parts of the universe. [In this case, the fluctuations 5(f)(1) on different 
scales play the role of the fields Xa] The question is then what spectrum 5p/p(l) we are 
most likely to observe in our neighborhood? In more general terms, we would like to deter- 
mine the probability distribution V(x) for the fields Xa associated with variable constants. 

The inflationary scenario implies a very large universe inhabited by numerous civilizations 
that will measure different values of Xa- We can define the probability V(x)dxi---dXn f° r 
Xa to be in the intervals dx a as being proportional to the number of civilizations which will 
measure Xa in that interval 0. This includes all present, past and future civilizations; in 
other words, it is the number of civilizations throughout the entire spacetime, rather than 
at a particular moment of time. Assuming that we are a typical civilization, we can expect 
to observe Xa near the maximum of V(x) [[|. The assumption of being typical has been 
called the "principle of mediocrity" in Ref. @. 

An immediate objection to this approach is that we are ignorant about the origin of 
life, let alone intelligence, and therefore the number of civilizations cannot be calculated. 
However, the approach can still be used to find the probability distribution for parameters 
which do not affect the physical processes involved in the evolution of life. The cosmological 
constant A, the density parameter Q and the amplitude of density fluctuations Q are exam- 
ples of such parameters. We shall assume that our fields Xa belong to this category. The 
probability for a civilization to evolve on a suitable planet is then independent of Xa, and 
instead of the number of civilizations we can use the number of habitable planets or, as a 
rough approximation, the number of galaxies. [We are assuming that civilizations can exist 
only for a finite period of time after inflation, either because the stars run out of nuclear fuel 
and die, or because protons eventually decay. A galaxy then gives rise to a finite average 
number of civilizations.] Thus, we can write 

V( X )d n x oc dAf, (1) 

where dhf is the number of galaxies that are going to be formed in regions where Xa take 
values in the intervals dx a - 

The meaning of Eq. (|IJ) is unambiguous in models where the total number of galaxies 
in the universe is finite. Otherwise, one has to introduce some cutoff and define the ra- 
tio of probabilities for the intervals d n x and d n x as the ratio of the galaxy numbers 
dj\f^ I dj\f^ in the limit when the cutoff is removed. The situation is relatively straight- 
forward in the case of an infinite universe which is more or less homogeneous on very large 
scales. One can then evaluate the ratio dN^'/dN^ 2 ' in a large comoving volume V and 
then take the limit as V — > oo. The result is expected to be independent of the limiting 
procedure; for example, it should not depend on the shape of the volume V. (It is assumed, 
however, that the volume selection is unbiased, that is, that the volume V is not carved to 
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favor some values of Xa at the expense of other values.) 

The problem is much more serious in models of eternal inflation where the universe 
is never completely thermalized ||[7j . In such models, the volumes of both inflating and 
thermalized regions grow exponentially with time and the number of galaxies grows without 
bound, even in a region of a finite comoving size. One can deal with this problem by 
introducing a time cutoff and including only regions that thermalized prior to some moment 
of time t c , with the limit t c — > oo at the end. One finds, however, that the resulting 
probability distributions are extremely sensitive to the choice of the time coordinate t [|J. 
This applies in particular to the calculations of the density fluctuation spectra. Having 
chosen t to be the proper time along the world lines of comoving observers, Linde, Linde 
and Mezhlumian || found that a typical observer could find herself near the center of a very 
deep minimum of the density field. On the other hand, if one uses the expansion factor along 
the trajectories as the time coordinate, one recovers the standard results §|. Coordinates 
in General Relativity are arbitrary labels, and such gauge dependence of the results casts 
doubt on any conclusion reached using this approach. 

An alternative procedure, suggested in |L(J, is to introduce a x-dependent cutoff at the 
time t c (x), when all but a small fraction e of the comoving volume destined to thermalize 
with Xa i n the intervals dx a has thermalized. The limit e — > is taken after calculating the 
probabilities. It was shown in [10,11 ] that the resulting probability distribution is essentially 
insensitive to the choice of time parametrization. However, the same problem appears in 
a different guise. Linde and Mezhlumian |12| have found a family of gauge-invariant cutoff 
procedures which includes the e-procedure described above. All these procedures give vastly 
different results for the probability distributions. 

In the face of these difficulties, serious doubts have been expressed that a meaningful 
definition of probabilities in an eternally inflating universe is even in principle possible 
PJ5|JT^]. Once again, this problem cannot be brushed aside as "exotica" because inflation 
is eternal in practically all models suggested so far. Eternal inflation is generic; it is finite 
inflation that is exotic. We believe, therefore, that the situation deserves to be called the 
"predictability crisis in inflationary cosmology". 

In this paper we are going to analyze the resolution of the crisis suggested by one of us 



in Ref. |13fl . We begin, in Section II, by discussing the spacetime structure of an eternally 
inflating universe. This will help us clarify the origin of the problem and will motivate its 
resolution proposed in Section III. In the following Sections IV - V we develop methods 
for calculating probabilities based on the prescription introduced in Section III. The first 
method uses the Fokker-Planck equation of stochastic inflation. In Section IV we derive 
the appropriate form of the equation, discuss the initial conditions, and work out some 
examples. The second method uses numerical simulations of eternal inflation. In Section 
V we describe our simulations and compare them to previous work. We also compare the 
results obtained from the simulations with those obtained using the Fokker-Planck equation. 
Our conclusions are summarized in Section VI. 



II. SPACETIME STRUCTURE AND GAUGE DEPENDENCE OF 

PROBABILITIES 

An inflating Universe can be locally described using the synchronous coordinates, 
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ds 2 = dr 2 — a 2 (x, r)dx 2 . 



(2) 



The lines of x = const in this metric are timelike geodesies corresponding to the worldlines 
of co-moving observers, and the coordinate system is well defined as long as the geodesies 
do not cross. This will start happening only after thermalization, when matter in some 
regions will start collapsing as a result of gravitational instability. Hence, the synchronous 
coordinates (0) can be extended to the future well into the thermalized region. 



A. Double-well model 

We shall first consider the simplest situation when there is a single inflaton field </> 
with a double-well potential, as in Fig. [I], and no other fields \a- The potential in the 
figure is symmetric with respect to <p — > —ip, although this is not essential. Depending on 
quantum fluctuations, an inflating region with ip « near the maximum of the potential 
can thermalize in either of the two minima at ip = ±77. The two minima may correspond to 
different physics, in which case this is an example of a model where the constants of Nature 
can take a discrete set of values. 

In the slow-roll regime of inflation, the field ip is a slowly-varying function of the coordi- 
nates, so that spatial gradients of <p can be neglected and 

y? 2 < 2V(cp). (3) 

The classical evolution of the scale factor a(x, r) and of the inflaton y?(x, r) is then described 
by the equations 

(a/a) 2 = H 2 ^ 87rV((p)/3, (4) 



« -V'(<p)/3H = -H'{y)/Air, (5) 

where dots represent derivatives with respect to r and we use the Planck units, h = c = 
G = 1. With the aid of Eqs. the slow-roll condition (|3|) can be expressed as 

H' < QH. (6) 

This condition is violated near the points (p = ip* \ signalling the end of inflation. The slow 
variation of <p implies that H is also a slowly-varying function of x and r, and thus the 
spacetime is locally close to de Sitter, with a horizon length if -1 . 

Quantum fluctuations of the field ip can be pictured as a 'random walk' [superimposed 
on the classical motion (JD] in which ip undergoes random steps of rms magnitude (5ip) rms = 
H/2tc per Hubble time, St = H~ l . The fluctuations Sip are correlated on the scale of 
the horizon (/ ~ H^ 1 ), but correlations rapidly decay with distance and become totally 
negligible on the scale of a few horizons. The fluctuations are dynamically unimportant if 
the classical 'velocity' \ip\ is much greater than the characteristic speed of the random walk, 
(^)rms/^ = H 2 /27T, which gives 

H'^H 2 . (7) 
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This condition is violated in the region tp)j-> < ip < (p@> near the top of the potential (see 
Fig. [I]). The dynamics of (p in this region is dominated by quantum fluctuations. The 
deterministic slow-roll regions are bounded by the values <p(p and (pi 3 , j = 1,2. 

The random-walk picture of quantum fluctuations can be used in numerical simulations 
of eternal inflation |H|,P| . To illustrate the spacetime structure of an eternally inflating uni- 
verse, we have performed simulations for several models in one and two spatial dimensions, 
with and without the additional fields Xa- The details of the simulations are given in Section 
V. Fig. [| shows the distribution of inflating and thermalized regions in two dimensions at 
four consecutive moments of time for the single-inflaton model with a double-well potential. 
The form of the potential is taken to be 

V{<p) = V cog i (a<p), (8) 

where a = 7i/2r]; we only consider the range —77 < <p < rj. With a ^ 1 and H Q /a < 1, we 
have 

^ ~ ±H /a\ (9) 
where H = H(<p = 0) = (87rV /3) 1/2 , and 

<pW « ±( V - 1/6). (10) 

In the simulations we used Vo = 0.05 and a = 1. Our choice of parameters in most of 
the simulations was dictated by the computational constraints (see Sec. V), so we made no 
attempt to make this choice realistic. 

There are two types of thermalized regions in the model ([5]) (corresponding to the two 
minima of the potential at (p = ±77) which are shown with different shades of grey. Ther- 
malized regions of different type can never merge; they are always separated by inflating 
regions (shown in white) which can be thought of as inflating domain walls ]T3| . This is an 
example of topological inflation. Thermalized regions of the same type can also be separated 
by inflating regions, but there is nothing to prevent these regions from merging, and we can 
see from the sequence in Fig. [| that neighboring regions of the same type do tend to merge. 

As time goes on, more and more comoving volume gets thermalized, and at r — > 00 only 
a vanishing fraction of the initial volume is still inflating. However, the physical volume 
of the inflating regions grows exponentially with time. Geometrically, these regions form a 



fractal of dimension d < 3 fll^ 



A spacetime slice through the universe (in the x — r "plane") is shown in Fig. |3|. The 
simulation starts with a horizon-size region and a homogeneous field (p = at r = (bottom 
of the figure). We shall be particularly interested in the thermalization hypersurfaces : 
(p = ip* which separate inflating and thermalized spacetime regions. We see from the figure 
that these surfaces tend to become asymptotically static in the comoving coordinates as 
t — > 00. This does not mean, however, that the surfaces become timelike at large r. To 
picture the geometry of these surfaces, one has to keep in mind that physical lengths in 
an inflating universe are related to coordinate differences by an exponentially growing scale 
factor. As a result, the boundaries of thermalized regions (as well as all surfaces of constant 
tp in the slow-roll regime) are spacelike, and are in fact very flat. 
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The spacelike character of the constant-^ surfaces can be understood as follows [10 
Consider the normal vector to the surfaces, d^ip. We have 



d^ip = (p 2 - a~ 2 (Vy?) 2 . 



(11) 



The spatial gradients of tp are caused by quantum fluctuations. On the scale of the horizon, 
the gradient is of the order a _1 |Vv?| ~ (S(p) Tms / H' 1 ~ H 2 /2tt, and is even smaller on larger 
scales. On the other hand, from Eqs. (§),(0), \<f \ ^> H 2 /2ti. Hence, d^ipd^ip > 0, and the 
surfaces <p> = const are spacelike. 

In general, thermalization hypersurfaces S* cannot terminate: they must either be in- 
finite or closed. A closed finite thermalization surface is possible only in a closed universe 
with finite (non-eternal) inflation. As Fig. [3] illustrates, thermalization surfaces in an eter- 
nally inflating universe extend indefinitely in the r direction and have, therefore, infinite 
volumes. It is easily understood that the number of disconnected thermalization surfaces in 
our double-well model is also infinite. Thermalized regions include inflating islands which in 
turn include thermalized regions of both types, and so on ad infinitum. Regions of the same 
type may later merge and be, therefore, bounded (in spacetime) by the same thermalization 
surface. But clearly there is going to be an infinite succession of nested regions of different 
types which are bounded by disconnected thermalization surfaces. 

For an observer in one of the thermalized spacetime regions, the surface S* at the bound- 
ary of that region plays the role of the big bang. The natural choice of the time coordinate 
in the vicinity of that surface is t = <p, so that the constant- £ surfaces are (nearly) surfaces 
of constant energy density. Since these surfaces are infinite, the observer finds herself in 
an infinite thermalized universe, which is causally disconnected from the other thermalized 
regions. The situation here is similar to that in the 'op en- universe' inflation, where ther- 
malized regions are located in the interiors of expanding bubbles and have the geometry 
approximating open (k = —1) Robertson- Walker universes. 

It is now easy to see why the probability distributions obtained using a cutoff at t — t c 
are so sensitive to the choice of the time variable t. Any spacelike surface S can be an 
equal-time surface t = t c with an appropriate choice of t. Depending on one's choice, the 
surface S may cross many thermalized regions of different types (e.g., for t = r), may cross 
only regions of one type, or may cross no thermalized regions at all (say, for t = ip with ip 
in the slow-roll range). These possibilities are illustrated in Fig. f| by surfaces Si,S 2 and 
S3, respectively. If, for example, one uses the surface <S 2 as the cutoff surface, one would 
conclude that all observers will see the same vacuum [the same minimum of V(<£>)] with 100% 
probability. With a suitable choice of the surface, one can get any result for the relative 
probability of the two minima. 



The potential is symmetric with respect to x ~ ¥ ~X> an d we shall not distinguish between 
the values x an d — X, so the effective range of the field xisO^X^ 71 "//^- The maximum of 



B. Two-field model 



We now consider a two-field model with a potential 



V(^p, x) = Vo cos 2 (aip)[l + A(l + cos fix) sin 4 atp}. 



(12) 
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the potential is at tp = and the true vacuum is at tp = ±77. The field x is niassless in the 
true vacuum, while the mass of tp is x-dependent, 

mj(x) = 2V a 2 [l + A(l + cos/?x)]- (13) 
The expansion rate in the model (|T^) depends on both tp and x> 

%x) = y%x)- (14) 

The quantum and thermalization boundaries, <p q and tp*, are also x-dependent; they are 
determined from 

W t .x)h%x), (15) 
I^Vv^x)! ~ 6iy(^»,x)- (16) 

An important special case is when H is essentially independent of % in the quantum regime 



and near thermalization. Then (p g and tp* are also nearly constant. Our model (|T2|) belongs 
to this class when A and (3 are sufficiently small. 

For small f3 and A, the x _ dependence of the potential is weak and the evolution of ip in 
the slow-roll regime is essentially the same as in the one- field model (|||), with \ remaining 
nearly constant. Different parts of the universe will thermalize with different values of % 
due to quantum fluctuations in the random-walk regime near the top of the potential. This 
will result in spatial variation of the amplitude of density fluctuations [|TJ, 

5p/p( X ) ~ 200t7v(x). (17) 

In the opposite regime, the classical evolution and quantum fluctuations of x are non- 
negligible during the slow-roll of ip. In this case, the density fluctuations will not be directly 
related to the value of x at thermalization, as in Eq. (|T7j). The quantitative conditions for 
the two regimes will be specified in Section IV. 

In order to make the quantum dispersion of x during the time of the simulation com- 
parable to the range of x-, we had to use a relatively large value of (3 = 10 (see Section 
V). We used the same values of Vo and a as before and A = 0.1. With these parameters, 
one finds that <p q and ip* are essentially independent of x and are approximately given by 
Eqs. ©,(0). 

The spacetime structure of the universe in this model is very similar to that in the one- 
field model (H). There is an infinite number of infinite thermalization surfaces, with the field 
X varying continuously along these surfaces. Fig. |5] gives a snapshot of a simulation where 
different values of x at thermalization are shown with different shades of grey, from light 
grey for x = to dark grey for x = 7r/ /3. Inflating regions are white as before, and here we 
do not distinguish between regions thermalizing at ip — +rj and ip = —77, assuming that they 
have identical physics. A spacetime slice through this simulation is shown in Fig. |6|, where 
now the values of x are indicated throughout the inflating region using the same shading 
code, while thermalized regions are white. 

Let us now try to understand the gauge dependence of probabilities obtained by a 
constant-t cutoff in this model. The simple explanation that we gave for the one-field 



7 



model, that constant-t surfaces can be chosen so that they cross one type of thermalization 
surfaces and avoid the other, does not apply here. The field x takes all its possible values 
on each thermalization surface E*, and the distributions of x on different surfaces should 
be statistically equivalent. However, when we change the time variable, say, from t = r to 
t = a, the shape of the cutoff surface t = t c is also changed, and it is important to note that 
the direction in which this surface is deformed is correlated with the local value of x- The 
scale-factor time a is greater in regions where the expansion rate H is higher, and since the 
expansion rate is x-dependent, the cutoff surface deformation will tend to favor regions with 
X corresponding to low values of H and will tend to exclude regions where x corresponds to 
high values of H. This is the origin of the gauge dependence of the cutoff procedure. 

One might have thought that the deformation of the cutoff surface due to the change of 
the time variable should affect only the boundaries of the thermalized regions included in the 
calculation of probabilities and should therefore have little effect in the limit t c — > oo when 
the volume of the thermalized regions gets very large. However, the volume of thermalized 
regions in an eternally inflating universe grows exponentially with time, and at any time 
the total thermalized volume is dominated by the newly thermalized regions which stopped 
inflating within the last few Hubble times.[] Hence, most of the thermalized volume is always 
in the vicinity of the cutoff surface t = t c , and it is not surprising that the probability 
distribution V(x) is sensitive to the choice of the time variable t used to implement the 
cutoff. 

The spacetime structure may be different for other models of eternal inflation, partic- 
ularly for topological inflation [15] with a different topology. For example, in the case of 
inflating monopoles one can expect inflating regions to be localized and to be surrounded by 
thermalized regions. It is conceivable that in this case there is a single connected thermalized 
region of spacetime which is separated from the inflating regions by a single thermalization 
surface £*. Despite this possible difference, we still expect that the constant-time cutoff 
procedure will still give gauge- dependent results, for the same reasons as in the two-field 
model discussed above. 



III. THE PROPOSAL 



We now review the resolution of the gauge dependence problem proposed in Ref. |T3| . 
Let us first assume that inflating and thermalized regions of spacetime are separated by a 
single thermalization surface £#. The problem with the constant-time cutoff procedures is 
that they cut the surface in a biased way, favoring certain values of x an d disfavoring 
other values. We thus need an unbiasly selected portion of £*. 

It appears that the simplest strategy is to use a "spherical" cutoff. Choose an arbitrary 



1 Rcmcmber that Figs. 5,6 show the comoving volume distributions of x- The physical volumes are 
related to the co-moving volumes by the factor a 3 (x, t), so the volumes thermalizing at later times 
are exponentially enhanced. Moreover, as Figs. 5,6 indicate, a fixed-time cutoff tends to produce 
a multitude of disconnected thermalized regions. Their total volume is dominated by small, newly 
formed regions. 
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point PonS,. Define a sphere of radius R to include all points Q whose geodesic distance^ 
from P along is d(Q, P) < R. We can use Eq. (HD to evaluate the probability distribution 
V(x) in a spherical volume of radius R c and then let R c — > oo. If the fields Xa var Y in 
a finite range, they will run through all of their values many times in a spherical volume 
of sufficiently large radius. We expect, therefore, that the distribution V(x) will rapidly 
converge as the cutoff radius R c is increased. We expect also that the resulting distribution 
will be independent of the choice of point P which serves as the center of the sphere. 

The same procedure can be used for fields with an infinite range of variation, provided 
that the probability distributions for Xa are concentrated within a finite range, with a 
negligible probability of finding Xa very far away from that range. 

Suppose now that there is an infinite number of disconnected thermalization surfaces, 
as in the two-field model of Section II. We can then pick an arbitrary connected component 
of £* and apply the spherical cutoff prescription described above. Since the inflationary 
dynamics of the fields Xa have a stochastic nature, the distributions of Xa on different 
connected components of £* should be statistically equivalent, and the resulting probability 
distribution V(x) should be the same for all components. 

An even more complicated spacetime structure is found in models with a non-vanishing 
cosmo logical constant. In such models, regions of true vacuum in the post-inflationary 
universe may fluctuate back to the quantum diffusion range of the inflaton potential. Every 
such region will serve as a seed for a new eternally inflating domain whose internal structure 
will resemble that shown in Fig. [5[ Thermalized regions formed in this domain will in turn 
produce new inflating seeds, etc. ||16|| . In such a "recycling" universe, there is an infinite 
number of thermalization surfaces to the future of any given thermalization surface. But once 
again, we expect them all to be statistically equivalent, so the prescription for probabilities 
remains the same. 

We note that this prescription cannot be applied to models where the potential has a 
discrete set of minima, as in the double-well model of Section II. We can introduce a discrete 
variable x labelling different minima. Each connected component of the thermalization 
surface E* is then characterized by a single value of Xi an d it is clear that the probability 
distribution for x cannot be determined by studying one such component .Q This may indicate 
that no meaningful probability distribution can be defined for a discrete variable in an 
eternally inflating universe. 

Having specified the cutoff procedure, the problem of calculating V(x) can now be split 
into two parts. The number of galaxies dAf(x) in Eq. (|I|) is proportional to the volume of 
the comoving regions where Xa take specified values and to the density of galaxies in those 
regions. The volumes and the densities can be evaluated at any time. Their product should 
be independent of the choice of this reference time, as long as we include both galaxies that 
were formed in the past and those that are going to be formed in the future. It will be 
convenient to evaluate the volumes and the densities at the time when inflation ends and 



2 If there is more than one geodesic connecting P and Q, d(Q, P) is defined as the smallest of all 
the geodesic distances. 

3 We assume that the discrete vacua cannot be connected by non-inflating domain walls. 
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vacuum energy thermalizes, that is, on the thermalization surface E*. Then we can write 

V{ X ) oc v(x)r*(x)- (is) 

Here, V*(x)d n x is proportional to the volume of thermalized regions where Xa take values 
in the intervals dx a , and v{x) is the number of galaxies that form per unit thermalized 
volume with cosmological parameters specified by the values of Xa- The calculation of v(x) 
is a standard astrophysical problem which is completely unrelated to the calculation of the 
volume factor V*(x)- Our focus in this paper will be on the volume factor. In the following 
sections we shall discuss some methods that can be used to calculate 



IV. THE FOKKER-PLANCK EQUATION 

A. Derivation of the Fokker-Planck equation 

A statistical description of fluctuating quantum fields in an eternally inflating universe 
can be given in terms of the probability distribution on equal time surfaces, P(p, Xi t)> which 
satisfies the Fokker-Planck equation. This description was suggested in Ref. ||, where the 
idea of eternal inflation was first introduced, and was further developed in Refs. |IT|-pO|,f4|,|2T 



Although this is an elegant formalism, previous attempts to use it for the calculation of 
probabilities gave ambiguous gauge-dependent answers P|,|2"2||7[] . Here we suggest a version 



of the Fokker-Planck formalism which allows an unambiguous calculation of 7-*(x) i n a wide 
class of models. 

The idea is that for p ^> p q , when quantum fluctuations of ip (but not necessarily of x) 
are negligible compared to its classical motion, the evolution of p at any given co-moving 
point is essentially monotonic in proper time and we can use ip as a new time variable, 
t = <p. Let us first assume that the x-dependence of the expansion rate H is insignificant near 
thermalization, so that p* is also nearly independent of x- Then the thermalization boundary 
p = <p* is a constant- "time" surface and the probability distribution V(x) is given simply 
by P(p*,x)- We shall now derive the Fokker-Planck equation for P(p,x)- That equation 
will be applicable in the regime p ^> <p q where the evolution of p is essentially classical (but 
without such restrictions on x), and will have to be complemented by appropriate boundary 
conditions (see Section IV, B). 

The classical drift velocity of x i n time t = p can be written as 



v 



dx X H\ 



x 



dp p H',, 



(19) 



where dots indicate derivatives with respect to the proper time r, H(p, x) is given by Eq. (|14| 
and H' x = dH/dx, etc. The quantum dispersion during the "time" increment dp is 



, H 3 , H 3 



(Sx q ) 2 = ^dr = ~^rdp = 2D(p, X )dp, (20) 



where D(p,x) is the diffusion coefficient. Finally, the rate of the cosmological expansion is 

Ida H A H 

-T- = - = -47T-— . 21 

a op p H' 
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The Fokker-Planck (FP) equation for P(<p, x) can now De written as 



OP d 2 d H 

^^ {DP) -a- x ^ p) - M w/- < 22 > 

where d is the number of spatial dimensions. The first and second term on the right-hand 
side of (^) describe the quantum diffusion and the classical drift of Xi respectively, and 
the last term accounts for the expansion of the universe. Eq. (B^) can also be derived 
directly from the FP equation for the two-field distribution function P(<p, x> t); this is done 
in Appendix A. We should note that due to the growth of physical volume, the distribution 
P(<p,x), just like P(<f,x,t), is not normalized. The physical probability distribution for x 
is given by P(<p*, x) after an appropriate normalization of that function over the range of x- 



Starobinsky has pointed out (17| that the diffusion term of the FP equation suffers from 
the ambiguity in the ordering of non-commuting factors d/dx and D(p,x)- For example, 
one could write this term as 



d_ 

dx 



£,1/2+7 A f D l/2-Vp 

dx \ 



(23) 



However, it has been recently shown |23| that the so-called Ito factor ordering, 7 = —1/2, 
has significant advantages compared to other choices. Hence, we used 7 = —1/2 in Eq. (p2|). 

It is also possible to relax the assumption that the thermalization surface = 0* is 
independent of x- Usually the endpoint 0* is defined only approximately and one could 
redefine its exact position (within the range where fluctuations of both x and are negligible) 
to remove the dependence of 0* on Xi without changing the physical results. In case such 
redefinition is impossible, the formalism can be extended to give the distribution of x along 
an arbitrary boundary line specified by a function 0* (x) ■ To implement this, one formally 
puts the diffusion coefficient D (0, x) and the drift velocity v x to zero everywhere in the 
domain > 0* (x) and solves the modified FP equation in the entire range of and x U P 
to a suitable value = max which should lie beyond the boundary 0* (x) for all x- The 
solution of the new FP equation adequately describes fluctuations of x up to the boundary, 
while beyond the boundary the distribution is kept unchanged for a given x- The restriction 
of this solution to = max (clearly independent of the choice of max ) gives the desired 
probability distribution V* (x). We shall therefore assume below that the boundary 0* is 
independent of X- 



B. Initial conditions 

To solve the FP equation ( p2|) we have to specify the initial distribution P(<po, x) at 
some <p = <p in the deterministic slow-roll range, ip q (x) <S < <P*(x)- We do not know 
how to do this in the general case. However, a simple initial distribution can be specified in 
models where the potential V(ip, x) is very symmetric near the top, so that it is essentially 
independent of x i n the interval < (p < <p\ which includes the diffusion range and part of 
the slow-roll range {ip q <C (fi < (p*). We can then choose <po to be in the part of the slow-roll 
range where the potential is still symmetric (<p q <C y^o < (pi). 

Let us consider the (spacelike) surface S : <p(x) = p . In the slow-roll regime <p rolls 
from smaller to larger values, and thus the values of <p in the past of S are in the range 
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<p q ?S <p < fo, and prior to that in the diffusion range — <p q ^ v 9 ~ W m both of these ranges 
the potential is nearly independent of Xi an d therefore we expect the distribution P((po, x) 
on E to be flat, 

P(<p , x) = const, (24) 

with all values of x equally probable. 

We shall now specify the conditions of applicability of the distribution (|24]) . Suppose the 
potential V(<p, x) can be represented as 

V(<p, X ) = U(<p) + V(<p,x) (25) 
with V< t/. The expansion rate can be written as 

H{<p,x) = H{<p) + 6H{<p, X ), (26) 

where 

«^, X )~47rV(^x)/3tf(^). (27) 

Let us consider the physical processes affecting the distribution of x an d their characteristic 
timescales. The x-dependent correction to the expansion rate (p7|) will tend to distort the 
distribution, producing a peak at the maximum of V(ip, x) where the expansion rate is 
the highest. This differential expansion becomes important on a timescale T de such that 
5Hr de ~ 1, 

r de ~ l/SH. (28) 

On the other hand, the classical force will tend to drive x to the minimum of V(<p, x) on a 
timescale 

r c ~^~47r^, (29) 
v c On 

where we have used v c = —5H' x /Att ~ 5H/A'kA x for the classical drift velocity.^ Finally, the 
quantum diffusion of x wm tend to keep the distribution flat. The typical variation of x due 
to the quantum random walk in a proper time interval r is (A%) 2 ~ (H 3 /4tt 2 )t. The time 
r q it takes x to spread over its entire range A x is 

r q - A 2 X /H 3 . (30) 

The effect of randomization due to diffusion will prevail over the differential expansion 
and classical drift, provided that 

Tq < T c , T de , (31) 



4 Note that v c is different from v x in Eq. (|8|) because the velocity v x is defined with respect to 
"time" tp, while v c is defined with respect to the proper time. 
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or, 



A15H 1 6H 

■in- h T.w « < 32 > 

If conditions (|32| ) are satisfied at tp < then we expect that P(<p g , x) ~ corasi at the 
onset of the slow roll. At <p 3> the inflaton ip grows monotonically with time, causing an 
evolution of the expansion rate H and of the potential V. This introduces two additional 
timescales 

r H ~ H/H ~ AitH/H'l (33) 

and 

-~ V/V ~*^V (34) 

If both of these times are large compared to r q , then the flat distribution of \ extends into 
the slow roll regime for as long as r q remains the smallest of all relevant times. In fact, 
the evolution of H and V does not by itself distort the flat distribution of \i an d one can 
expect that this distribution persists as long as the conditions ([32|) are satisfied, even if r q 
gets larger than t# or ry. 

Another interesting case where the probability distribution for x can be specified is when 
r c and r q are small compared to all other relevant timescales, 

T c , T q < T de , T H , T V . (35) 

Then the distribution P(<p, x) is determined by a statistical equilibrium between the quan- 
tum diffusion and the classical drift. This equilibrium distribution is given by P|| 



P eq (<p, x) oc exp I 3 ff4(y,) J • ( 36 ) 



It is an approximate solution of Eq. fl22|) when the (^-dependence of D and v x can be disre- 
garded and the last term in (f22|) can be dropped. We note that it follows from Eqs. (P^ , (|3T)|) 
and (|3) that 

87r2y ***** (vr\ 

-3H*~ 87l 7 c - (37) 

Hence, for 8n 2 T q <C r c the distribution (^) is nearly flat, as expected. In the opposite limit, 
t~c <C r q , the distribution is strongly peaked at the minimum of V(<p, x)- If the conditions of 
equilibrium (^) are satisfied at some y?o in the range <C <po < we can use P eq (ipo, x) 
as the initial condition for the FP equation at (p — tpo- 
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C. Some examples 



We shall now give some examples of numerical and analytic solutions of the FP equa- 
tion ( p^ ) for the two-field model (|T2|). In this model U(<p) = V cos 2 (a(p), V((p,x) — 
XV cos 2 (aip) sin (aip)(l + cos/3x), A x = n/f3, and the conditions (i32|) which ensure a flat 
distribution at the onset of the slow roll take the form 

\H 2 /8na 4 < 1, (n/(3) 2 (\H 2 /a 4 ) < 1. (38) 

Here we have used (p ~ (p q ~ H /a 2 . 
In our first example we used 

a = 0.1, 13 = 0.1, A = 0.01, H = 10- 6 (39) 

and the initial condition P = const at <^ = 10 -3 (which is above ip q ~ 10~ 4 ). These 
values of the parameters could occur in a realistic inflationary model: Eqs. (|13"D and (|17D 
give ~ 5 x 10 -7 and Sp/p ~ 10~ 5 . The characteristic times for the parameters (|39| ) are 
plotted in Fig. |7| as functions of ip. We see that the conditions (|31~D are initially satisfied, 
and thus we are justified to use a flat distribution as the initial condition. The result of a 
numerical integration of Eq. (p2|) is shown in Fig. [8]. At <p> = <p>*, the distribution is peaked 
at x — 0, indicating that the effect of differential expansion was more important than that 
of classical drift and of diffusion. This is not surprising, since r^ e C r,, r c in most of the 
slow roll range. 

The classical variation of \ during the slow roll of tp from some ipo > f q to ip* is given 

by 

$Xc = r ^rdtp ~ X(3/a 2 (40) 

and is essentially independent of <p> . Similarly, the variation of x due to quantum fluctuations 
is 

, x n 2 f H 3 , 1 r* R3 a H o i sin(a^) 

J A.tt z 7r Jipo H'lp not 1 sm(aip ) 

with only a logarithmic dependence on <po. For the parameter values in our example, both 
5xc and 5x q are small compared to the range of Xi = 71 /P, and thus x remains nearly 
constant during the slow roll. 

When both diffusion and classical drift are negligible, the first two terms on the right- 
hand side of the FP equation fl22|) can be dropped, and the solution of ( p2|) is 

P(,,,,)«exp(-«/;|I^,,)^^. (42) 

Here, 

N( X ) = I ' Hdr = -4vr f * JLd<p ^N + ^ cos(/3 X ) (43) 
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is the number of inflationary e-foldings between <po and </?*, Nq is the number of e-foldings 
in the ^-independent case (A = 0), and we have used sin(ay9^) « 1, which is valid for a ^$ 1 
[see Eq. QTO)] . Eq. (|42| ) has a simple interpretation [[13]: if x remains nearly constant, then 
the probability distributions for x at V 9 = fo an d V 9 = V 9 * are related simply by the volume 
expansion factor e . Our numerical solution (which was done for d = 3) agrees with the 
analytic solution fl42|) within 10%. 
Our second example is 

a = 0.01, /3 = 10 4 , A = 10~ 8 , H = 10" 4 (44) 

with a flat initial condition at tp = 10. The characteristic timescales for this case are plotted 



in Fig. [5] and the numerical solution of the equation is shown in Fig. |TU[ In this example, r q is 
initially much smaller than all other characteristic times, and thus the use of the flat initial 
condition is justified. The distribution develops a peak at % = ir/j3 at "time" <p ~ 60, when 
87r 2 Tg/r c ~ 1, as expected [see Eq.([37])]. r c and r q remain smaller than other characteristic 
times throughout the slow roll period, so one expects that the equilibrium distribution ( |36|) 
should be an approximate solution of the FP equation (|22|). We have verified that this is 
indeed the case, with a very good accuracy. The two distributions begin to diverge as we 
approach thermalization (<p ~ 140), presumably because ry becomes comparable to r c and 
r q . However, even at <p = 140 the equilibrium distribution is accurate within ~ 5%. We had 
to terminate the numerical solution at tp = 140, somewhat short of <p* ~ 157, because the 
peak of the distribution was getting too narrow. 

Finally, we present a numerical solution of Eq. ( |2"2| ) for the parameters we used in the 
simulations: 

a = 1, p = 10, A = 0.1, H = 0.63. (45) 

with a flat initial condition at p>o = 0.8. The simulations were performed in (1 + 1) di- 
mensions, so we used d — 1. (This choice of parameters was dictated by computational 
constraints; see Section V). The characteristic timescales for this example are plotted in 



Fig. [TT| and the solution is shown in Fig. [12]. Once again, we see that initially r q is smaller 
than all other timescales, and thus our choice of the initial condition is justified. But to- 
wards the end of the slow roll, r c , r de and t h become comparable to r q , within an order 
of magnitude. So in this example the shape of the final distribution at (p — (p* is hard to 
predict without a numerical solution of the FP equation. 

For the parameter values fl45|) we have <p q ~ 0.63, which is close to the initial value 
(po = 0.8. Strictly speaking, (p cannot be used as a time variable in such proximity of 
ip q . However, our solution of the FP equation is in a good agreement with the distribution 
P((p, x) obtained from numerical simulations (see Section V). This indicates that our method 
is robust and gives a reasonable accuracy even for p>o ~ <p q . 



V. NUMERICAL SIMULATIONS 

The most direct approach to the calculation of V*(x) is to use a numerical simulation 
of eternally inflating spacetime. The first simulation of this sort JL4| was performed on a 
square lattice representing a comoving region of initial size H~ l in a (2 + l)-dimensional 



15 



universe. The initial value of the inflaton ip was set to <p = at all points of the lattice. 
The potential was assumed to be very flat in the whole range \tp\ < tp*, so that H « const 
throughout the inflating region. Then there is no classical slow roll, and the evolution of 
ip is driven entirely by quantum fluctuations. The fluctuations were simulated by adding 
random increments to (p independently in each horizon-size region. After each two-folding of 
expansion, square regions which had physical size H^ 1 at the previous step were subdivided 
into four smaller squares and a fluctuation, drawn from a Gaussian distribution of width 
(if/27r)(ln2) 1//2 , was added to <p in each of these squares. Whenever \ip\ exceeded the 
corresponding cell was marked as thermalized. 

Linde, Linde and Mezhlumian f25|j4}j developed more realistic simulations which allow 
for the slow roll of <p. They performed simulations of inflating domain walls and vortices 
in (2 + 1) dimensions. The classical change of ip in a small proper time increment At was 
taken as [see Eq. @] 

Aif = ~^w^ At (46) 

and quantum fluctuations were represented by sinusoidal waves 

TT 

<fy,(x) = -^VHAt sin (He Hr {n ■ x + a)) (47) 



V27T 

of wavelength comparable to the horizon, I = 2n/ H. A phase a and a unit two-dimensional 
vector n were chosen at random at each step of the simulation. The expansion rate H was 
still treated as a constant. 

LLM have also performed (1 + 1) and (2 + l)-dimensional simulations in which they 
included the effect of variation of the expansion rate H in space and time [|J. In (1 + 1) 
dimensions they still used the sinusoidal ansatz for the fluctuations, but now the amplitude 
and the wavelength of the sine depended on the local value of H(x, r). In (2 + 1) dimensions 
with a variable H(x, r) this ansatz could not be used and LLM represented the fluctuations 
by wavelets. 

An advantage of using the sine waves is that they give a continuous variation of <p, in 
contrast to the discontinuous jumps along the boundaries of the squares in the method 



of Ref. | I4]| . However, use of the sine functions presents two problems. Firstly, the sine 
waves introduce unphysical correlations between the values of 5<p at points separated by 
large distances^]. In this respect the simulation of Ref. [14], which treated 5<p in neighboring 



horizon regions as completely uncorrelated, was more realistic. We note, however, that the 
calculation of probability distribution will be unaffected by these extra correlations, provided 
that the values of x are sufficiently well sampled throughout the volume. Also, the wavelet 
representation of fluctuations, as used by LLM in 2 + 1-dimensional simulations, eliminates 
correlations at distances beyond the wavelet size (which was taken to be larger than the 
horizon size in LLM). 



5 The unphysically large long-distance correlations resulting from Eq. ([17]) disappear only in the 
limit of At — > 0. These issues will be addressed in more detail in Ref. Il2l 
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Secondly, the fluctuation amplitude computed with the sine wave ansatz is strictly 
bounded from above by the value it takes when the sine function in Eq. (|47D is equal to one. 
In reality, the amplitude of the fluctuations is a Gaussian variable, and large values of 5tp are 
suppressed but not strictly forbidden. LLM have dealt with this problem by using a small 
timestep (At <C H~ x ), effectively using a superposition of many sine waves with random 
phases in each Hubble time. This approximates the desired Gaussian distribution signifi- 
cantly better than a single sine wave (although it does not eliminate unphysical correlations 
at large distances). 

Motivated by these considerations, we have attempted to develop simulations that do 
not suffer from these problems, as described in Sections V A. 

Numerical simulations of eternal inflation encounter significant computational con- 
straints, and we were able to perform simulations only for a restricted range of parameter 
values. Suppose we do the simulation starting with ip = in a horizon-size region, / = Hq 1 . 
After N e-foldings of expansion, the size of our comoving region is / = H ~ 1 e N . The simula- 
tion has to stop when there are still a few grid points per horizon, and realistically the total 
number of points cannot exceed M m ax ~ 10 7 . Then, in a (d + l)-dimensional simulation, the 
number of inflationary e-foldings is 

N < dT 1 \nAf max + ln(iV#*) ~ 16/d, (48) 

where we have assumed that ln(Ho/H*) ~ ln(6/a) ^ 16/d. 

The number of inflationary e-foldings in the slow-roll range (f q < if < tp* is 

.^fl^^^hf (49) 
Jtp q H \p cr sm(a(p q ) cr Hq 

and it follows from Eq. ([S|) that we need a ^ 1. On the other hand, a is bounded from 
above by a ^ 1, since otherwise there is no eternal inflation in this model. We used the 
value a = 1. 

The dependence of N on Hq is only logarithmic, but still the number of e-foldings is too 
large for small values of Hq. We used H = 0.63 (which corresponds to V$ = 0.05). This 
value of H is dangerously close to the Planck scale, but here we disregard all complications 
of quantum gravity. (Our choice of parameters in these simulations was dictated by the 
computational constraints, and we made no attempt to make this choice realistic.) 

In two-field simulations, the dispersion of x due to quantum fluctuations is 

5 Xq ~ « O.Q^= (50) 

^ Vd 

In order to have a representative simulation, we require that x should spread through its 
entire range (vr//3) by the end of the simulation. This implies 5x q ~ 7r//3, or (3 ^ 5d l ' 2 H~ 1 . 
In our simulations we used d — 1, /3 = 10. 

A. Comoving space simulations 

The first series of simulations described the spacetime in co-moving coordinates. The 
simulation was performed on a co-moving grid of points in either 2 or 1-dimensional space. 
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Both the one-field and the two-field models were simulated. Initially the field values were 
taken constant throughout the grid, with ip = 0. At each timestep, the fields were incre- 
mented by the classical evolution term Aip from Eq. ( ^6] ) and by the corresponding quantum 
fluctuation 5(p (x) . The evolution of the field x m the two-field model was computed in the 
same manner. The simulation was continued until the physical distance between adjacent 
grid points grew larger than the horizon size.[] The expansion factor was calculated for each 
point to determine the probability distribution of \ weighted by the physical volume. 

Instead of taking the sine waves, we used for 5(p (x) true Gaussian random fields with 
an appropriate correlation function. To obtain the increments for the grid points at each 
step, one only needs to generate a set of Gaussian variables dtp (xj) with zero mean and 
correlation matrix 



(Sip (x^ 5<p (xj)) = Qj = C (d (x^ x,-)) , (51) 

where d is the geodesic distance between points x«, x^ calculated using the known scale 
factors at intermediate points. Given the symmetric and positive-definite correlation matrix 
Cij, one can compute its Cholesky decomposition C = LL T where L is a lower triangu- 



lar matrix (see e.g. Ref. [|S7| for a numerical method), and the required vector of random 
increments is then obtained as 

5^ (Xi) = L\uj (52) 

where u is a vector of uncorrelated random numbers with normal distribution. The correla- 
tion function C(r) decays rapidly over the distance of a few horizon sizes and was taken to 
be zero at distances larger than 2H~ l to simplify calculations. 

We have performed simulations for the double-well model (^) and the two-field model 
(O). Momentary snapshots and spacetime slices of some of these simulations are shown in 



Section II. A more systematic investigation of the spacetime structure of eternally inflating 
universes, based on numerical simulations of this sort, will be given in a separate publication. 



B. Physical space simulations 



In the second series of simulations we have investigated the thermalization surface of 
the two-field model of Eq. ( |T2"D with the parameter values ( [45] ) in (1 + 1) dimensions. We 
have already discussed the spacetime structure of the universe in this model (Section IIB). 
As illustrated in Fig. |5], the universe contains an infinite number of disconnected, infinite, 
spacelike thermalization lines (thermalization surfaces are lines in a (1 + l)-dimensional 
simulation). In order to implement our prescription for the calculation of probabilities, we 
have to choose one of these lines and calculate T*(x) d£/dx on a large segment of that 



6 In the (2 + l)-dimensional case, we extended the simulations somewhat beyond this point and 
allowed the separation of adjacent grid points to exceed the horizon in the inflating regions with 
<p near the top of the potential. This should not have had any effect on the thermalized regions 
shown in Figures 2 and 5. 



18 



line. Here, d£ is the length of all parts of the segment where \ takes values in the interval 
dx- A segment of length 2R is a one- dimensional analogue of a sphere of radius R. The 
probability distribution is obtained in the limit of large R, so we used the whole length of 
the thermalization line that we were able to generate. 

The simulation space is a line of gridpoints, each carrying a value of <p and \\ initially the 
values of the fields are set constant everywhere. At every timestep, the fields are evolved by 
adding the classical drift Acp and the quantum fluctuation 5ip (x) generated from Eq. (E7p 
(we chose this over the more precise Eq. (|5^ ) for reasons of computational feasibility) in 
which the sine wave was multiplied by a random factor M, 

TT 

6<p(x) = M^^VHAt sin (He Hr {n ■ x + a)) . (53) 

It can be verified that if a is uniform in [0, 27r] and M has the standard x 2 distribution 
with two degrees of freedom, then Msma has normal distribution. This modification of 
Eq. (|47|) circumvents the problem we mentioned above with the unphysical distribution of 
the fluctuation amplitude resulting from a fixed-amplitude sine wave. 

The local expansion factor is also computed at each point. Due to the one-dimensional 
geometry of the simulation space, it is possible to maintain an (approximately) constant 
number of points per physical horizon size by inserting new grid points as necessary at each 
step. Thus we obtain a connected piece of the thermalization line in physical (instead of 
co-moving) coordinates. 

Since the calculation involves only a single thermalization line, one can save a great deal 
of memory and computer time by discarding parts of the simulation that do not affect that 
particular line. At any given moment the universe in the simulation consists of thermalized 
regions of two types (corresponding to tp = ±77) separated by inflating regions. Suppose 
we have chosen a plus-region (with ip = +rf) for the calculation of V*(x)- This region will 
never merge with any of the minus-regions, so we can safely discard the nearest minus- 
regions on both sides and everything beyond them. New minus-regions can later be formed 
in the remaining part of the simulation, and we can discard them and everything beyond 
them again. A greater efficiency is achieved by discarding the nearest thermalized regions 
regardless of their type, so the simulation includes only one thermalized region and two 
adjacent inflating regions. Once again, the procedure is repeated whenever new thermalized 
regions are formed. If the discarded nearest neighbors were plus-regions, they can later 
merge with our region of choice, in which case the simulation will terminate prematurely. 
This can be fixed by going back to the moment when the merged neighboring region was 
discarded and performing the simulation for that region from that time on. The resulting 
thermalization line is then attached to the line obtained for the initial region. We can make 
an arbitrary number of such attachments to generate a thermalization line of any desired 
length. We used this attachment procedure in our simulations. 



Fig. [L3| shows the probability distributions we obtained for two different thermalization 
lines in our simulations. As expected, the two distributions are very similar. Also shown 
is the distribution we obtained by solving the Fokker-Planck equation in Section IV C. It 
agrees with the simulation-based distributions with an accuracy of 4% . 
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VI. CONCLUSIONS 



In this paper we analyzed the origin of the gauge dependence of probability distributions 
for observables in models of eternal inflation. We reviewed the resolution of this problem 



proposed in Ref . |L3[ , gave a more precise formulation of this proposal and used it to develop 
two methods for calculating probabilities. 

The root of the problem is that an eternally inflating universe contains an infinite number 
of observers (galaxies), and one needs to introduce some cutoff procedure in order to calculate 
probabilities. The proposed procedure can be summarized as follows. Pick an arbitrary 
galaxy Q at an arbitrary moment in its history; this defines a point in spacetime. Consider 
the hypersurface of constant energy density S passing through that point. Define a sphere 
of radius R to include all points whose geodesic distance from Q along £ is less than R. 
The proposal is to evaluate the probability distribution P(x) based on galaxies formed 
within the co- moving region defined by a sphere of radius R c and then let R c —>■ oo. The 
surface £ is an infinite spacelike surface, and the distribution of the fields Xa on this surface 
is statistically homogeneous on large scales. One expects, therefore, that the resulting 
probability distribution P(x) will be independent of the choice of the "central" galaxy Q and 
that it will rapidly converge once the cutoff radius R c becomes larger than the characteristic 
scale £ of variation of the fields Xa on E. 

One argument in favor of this proposal is that it satisfies the "correspondence principle" 
with models of finite inflation, where the definition of probabilities is unambiguous. In such 
models, the total number of galaxies in the universe is finite, and all galaxies should in 
principle be included in the calculation of P{x)- However, a spherical cutoff with R c ^> £ 
should give a distribution very close to that obtained using the complete set of galaxies. 

We performed numerical simulations of eternal inflation in a two-field model ([12]) and 
showed how the spherical cutoff can be implemented to calculate the probability distribution 
P(x) numerically. Due to significant computational constraints, we were able to perform 
simulations only for a rather restricted set of parameter values. 

An alternative approach, which does not suffer from these restrictions, is to calculate 
P(x) by solving the Fokker-Planck equation with the inflaton field p playing the role of 
time. This form of the FP equation is valid only in the slow-roll range, ip q < p < ip*, and 
the method can be used only in cases where one can specify the initial distribution P(<po, x) 
at some po ^ (p g . We have shown in Section IV that this can be done in a reasonably wide 
class of models. In more general models, one may have to do numerical simulations in the 
quantum diffusion regime in order to determine the initial conditions for the FP equation 
in the slow-roll regime. 

We have indicated some cases in which the FP equation can be approximately solved 
analytically. In one of these cases, when the fields Xa remain nearly constant during the 
slow roll of p, the resulting distribution agrees with that obtained in We obtained 

a numerical solution of the FP equation for the parameter values used in our simulations 
and verified that the probability distribution P(x) obtained in this way agrees with the 
distribution calculated directly from the simulations. 

The approach we developed in this paper can be extended to models of "chaotic" inflation, 
where the inflaton potential rises to super-Planckian values. However, in such models one 
has to deal with the uncertainties associated with Planck physics. In particular, one has to 
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impose some boundary condition for the FP equation at the "Planck boundary" ip = tp p , 
where V(ip p ) ~ 1. 

Perhaps the most important application of our approach would be a probabilistic analysis 
of the density perturbation spectra. As we emphasized in the Introduction, this problem 
arises in practically all inflationary models, with or without the additional fields \ a - A 



brief discussion of this issue can be found in Ref. |T3[ where it was argued that the present 
approach should give essentially the same results as the standard calculations ||. A more 
detailed analysis will be given elsewhere. 

Our prescription for calculating probabilities cannot be applied to models where \ is 
a discrete variable. In such models, different values of \ are taken in different, causally 
disconnected thermalized regions. The spherical cutoff prescription involves only a single 
thermalized region, and it is clear that it cannot be used to determine the probability 
distribution in this case. One can take this as indicating that no probability distribution for a 
discrete variable can be meaningfully defined in an eternally inflating universe. Alternatively, 
one could try to introduce some other cutoff prescription to be applied specifically in the 



case of a discrete variable. It is possible that some version of the e-prescription of Ref. (TC 
could be used for this purpose. This issue requires further investigation. 
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APPENDIX A: DERIVATION OF THE REDUCED FOKKER-PLANCK 
EQUATION IN THE LIMIT OF NEGLIGIBLE DIFFUSION 

We start from the usual time-dependent FP equation in two dimensions (0, x) for the 
co-moving probability distribution P (t, 0, x), under assumption that diffusion in direction 
is negligible and using the Ito factor ordering: 

d t P = d\ (p (0, X ) P) - d x (v x P) - d+ {v+P) . (Al) 

Here and v x are the usual (both and x-dependent) drift coefficients, 

V4> = ~^ d 4> H > *x = ~^ d x H , ( A2 ) 

and the diffusion coefficient is D = H 3 /8tt 2 . The diffusion is bounded by the line = 0* 
which is the absorbing boundary (due to thermalization). No other exit boundaries are 
present. 

First we consider the probability density P e (t; 0, x\ 4>*i X*) to exit at a particular time 
t through the neighborhood of a particular point x* along the line = 0*, given that we 
started at t = at a given point (0,x)- Obtaining the distribution P e (t; ...) is a standard 
problem in the theory of stochastic processes [^HJ; once we know P e (t; ...), we can integrate 
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it over time until t = oo and obtain the total probability to thermalize in the neighborhood 
of X*- 

According to the standard theory, we can find the generating function g(s; ...) for the 
distribution P e (t; ...), defined as 



oo 

St 



g{s;<f>,x;<f>*,x*) = / e s P e (t;<f>,x;<f>*,x*)dt ( A3 ) 

Jo 

because it turns out that g is the solution of the "backward" or adjoint diffusion equation 

D (0, x) d 2 x g + v x d x g + v+d+g = sg (A4) 
with boundary condition at = 0* given by 

gis-^^x'A^x*) = Hx- x*) ■ (A5) 

This boundary condition makes the function vanish everywhere on the exit boundary except 
at the point x* where we want to find the exit probability. It is seen from Eq. (|A3|) that the 
total probability of exit at any time, i.e. the integral of P e through all t, is equal to the value 
of g at s — 0. This value g (s = 0; 0, x\ <f>*-> X*) = 9o is then interpreted as the probability 
that the diffusion process finishes at x — X* if if started at (4>,x)- The boundary condition 
Eq. ( |A5|) means that if we already started at the boundary but at a different Xi then the 
probability to exit at the chosen value x — X* is zero. The equation for g is 

D (0> X) 9 x g + v x d x9o = -Vfdfgo. (A6) 

This is formally the same as the backward Fokker-Planck equation for a one-dimensional 
diffusion process with <fi playing the role of time, where g (0, x] 4>*i X*) is interpreted as 
the probability of finishing at x = X* & t t ne final "time" <fi = 0* if started at (<f>,x)- The 
boundary conditions (|A5|) correspond exactly to that interpretation. The drift and diffusion 
coefficients become 

v x = ~% D = -° (A7) 

Vcj, V<f> 

in agreement with Eqs. (p!9|)-(p0|). The equation adjoint to Eq. ( |A6| ) is therefore interpreted 
as the diffusion equation for the "time" -dependent probability distribution P (0, x', Xo), with 
as "time" and the new kinetic coefficients, 

d^P = d\ (DP) - d x (v x P) . (A8) 

The equations for the physical volume distribution are derived in the same manner 
and differ from Eqs. ([Alp - flABD only by the addition of the growth term 3HP or 3Hg 
as appropriate, and the corresponding growth term in Eq. (|A8|) is similarly obtained as 
—SHP/Vfj). We therefore arrive at Eq. fl2"2]). 
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FIGURES 







V 

\ 
































W 1 ' 





FIG. 1. The inflaton potential for the double-well model (||). 




(c) (d) 



FIG. 2. A two-dimensional simulation for the double-well model at four consecutive moments 
of proper time: t = 5Hq 1 (&), t = 5.5-ffJ~ 1 (b), t = 6ff " 1 (c), t = 6.5ff ( ^ 1 (d). We evolved a comoving 
region of initial size I = Hq 1 with the initial value of f = at t = 0. Inflating regions are shown 
white, while thermalized regions with tp = +rj and ip = —rj are shown with different shades of grey. 
Thermalized regions of the same type tend to join in the course of the simulation. For example, 
regions labeled A and B in snapshot (c) have merged into a single region in snapshot (d). 
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FIG. 3. Spacetime structure in a one-dimensional simulation for the double-well model. It can 
be thought of as a spacetime slice through the (2 + l)-dimensional simulation illustrated in Fig. 2. 
Inflating regions are white, and thermalized regions of different type are shown with different shades 
of grey. 




x 

FIG. 4. S\ is a surface of a constant proper time. It crosses many thermalized regions of 
different types. S2 is a spacelike surface which crosses regions of only one type. 1S3 is a spacelike 
surface which does not cross any thermalized regions. 
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FIG. 5. A snapshot of a simulation for the two-field model (12). Inflating regions are white, 
while regions that thermalized with different values of x are shown with different shades of grey. 
The shading code is indicated in the bar on the right of the figure. 




X 

FIG. 6. Spacetime structure in a one-dimensional simulation for the two-field model showing 
the evolution of a comoving region of initial size I = Hq 1 with initially homogeneous (p = and 
X = 7r / /?. The values of the field x are indicated throughout the inflating region using the same 
shading code as in Fig. [|. The thermalized regions are left white. 
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FIG. 7. The characteristic times as functions of (p for the model parameters 
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FIG. 13. The probability distribution for \ f° r the model parameters (45). The distributions 
obtained directly from two different thermalization lines in a simulation are shown by dotted and 
dashed lines. The solid line shows the distribution obtained by solving the FP equation. 
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